//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "statics.hpp"
#include "screen.hpp"
#include "dynamics.hpp"

// test for error of b field in the center and along axis of a circular current
// using finite element analysis with different discretion sizes

#define RADIUS 1.0
#define CURRENT 5000
#define DIST1 0.1
#define DIST2 0.65
#define DIST3 3.65
#define DIST4 8.65
#define B_CENTER 0.0031415926535897932384626433832795028841971693993751058209747500
#define B_DIST1 0.0030950510166458728738817456004112390820202559707235045453751641
#define B_DIST2 0.0018517035538976034801836234486435147836112774342817613518642229
#define B_DIST3 0.0000579590991015778697875928556046998830675217920940527693531920
#define B_DIST4 0.0000047583098957210537076402953324873386690216600975446008312888
#define ALLOWED_FACTOR 1000.0 // allow difference by 1/1000 worth

int main (int argc, char *argv[])
{
  prec_t values[][2] = {
      {0, B_CENTER},
      {DIST1, B_DIST1},
      {DIST2, B_DIST2},
      {DIST3, B_DIST3},
      {DIST4, B_DIST4},
      {0, 0}
  };

  prec_t acc = 0.1;
  while ( acc > 0.00005 )
    {
      Statics::set_ddsq (acc);
      Statics statics;
      ElectroDynamics dyn(statics);
      statics.addStaticElement (new StaticRingIntegrated(5000, vect3d(0, 0, 0), RADIUS, 0, vect3d(0.0, 0.0, 1.0)));

      for ( size_t i = 0; values[i][1] != 0.0; i++ )
        {
          vect3d test(0, 0, values[i][1]);
          vect3d bf;
          dyn.getBField (vect3d(0, 0, values[i][0]), bf);
          vect3d vahe = bf - test;
          char b[128];
          prec_t perc = 100 * vahe.length() / values[i][1];
          printf ("ACC=%.7f, d=%.2f, err=%g (%.e%%), b=%s / (0, 0, %g)\n",
              prec2double(acc), prec2double(values[i][0]), prec2double(vahe.length()), prec2double(perc), bf.sprint (b), prec2double(values[i][1]));
        }
      acc = acc * 0.5;
    }

  return 0;
}

